clear
load('counter.mat')

theta_vector = [theta-0.25,theta-0.125,theta,theta+0.125,theta+0.25];
moments_vector = zeros(1,length(theta_vector));

for tthh = 1:length(theta_vector)
    
    theta_aux = theta_vector(tthh)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% STEP 1: Initial period (BGP) %%%%%%%%%%%%%%%%%%

data_pat_aux = data_pat(:,:,first_period) + data_pat(:,:,first_period-1) + data_pat(:,:,first_period-2);

u_0 = data_pop(:,first_period).^(1/zeta);
u_0 = u_0 / geometric_mean(u_0);
L_y_0 = (1/3)*(1/S)*repmat(data_pop(:,first_period),1,S);

lambda_0 = data_pat_aux ./ L_y_0;
lambda_0 = lambda_0 / mean(mean(lambda_0));

error_0 = 1000;

while error_0 > 0.0001
   
    kernel_utility_0 = u_0.*(lambda_0.^(1/theta_aux)); % The experience premium should appear as A_o^(1-betta) but does not matter, since it appears at the denominator too
    
    pi_0 = zeros(N,N,S);
    
    for j=1:N
        denom = sum(sum((kernel_utility_0./repmat(mu(j,:)',1,S)).^zeta));
        % Probability of going from j to any destination (:):
        pi_0(j,:,:) = ((kernel_utility_0./repmat(mu(j,:)',1,S)).^zeta)./denom;
    end

    error_L_y_0 = 1000;
    while error_L_y_0 > 0.1
       L_y_0_aux = zeros(N,S);
       for j=1:N
           L_y_0_aux = L_y_0_aux + sum(L_y_0(j,:))*squeeze(pi_0(j,:,:));
       end
       error_L_y_0 = max(max(abs(L_y_0 - L_y_0_aux)));
       L_y_0 = L_y_0_aux;
    end
    L_o_0 = L_y_0; L_k_0 = sum(L_y_0,2);
    pop_0 = L_k_0 + sum(L_y_0,2) + sum(L_o_0,2);
    
    lambda_1 = zeros(N,S);
    for i=1:N
        for s=1:S
            if L_y_0(i,s) > 0 && data_pat_aux(i,s) > 0
                lambda_1(i,s) = data_pat_aux(i,s) / L_y_0(i,s);
            end
        end
    end 
    lambda_1 = lambda_1 / mean(mean(lambda_1));
    
    u_1 = ((data_pop(:,first_period)./pop_0).^(1/zeta)).*u_0;
    u_1 = u_1./geometric_mean(u_1);
    
    error_0 = max(abs(u_0 - u_1)) + max(max(abs(lambda_0 - lambda_1)));
    
    rrr = rand;
    if rrr < 0.01
        error_0
    end
    
    lambda_0 = 0.99*lambda_0 + 0.01*lambda_1;
    u_0 = 0.99*u_0 + 0.01*u_1;
    
end
L_k_dyn(:,first_period) = L_k_0; L_y_dyn(:,:,first_period) = L_y_0; L_o_dyn(:,:,first_period) = L_o_0;
income_py_dyn(:,first_period) = sum(L_y_0.*(lambda_dyn(:,:,first_period).^(1/theta_aux)),2) ./ sum(L_y_0,2);
lambda_dyn(:,:,first_period) = lambda_0;
aggr_income_py_dyn(first_period) = sum(sum(L_y_0.*(lambda_dyn(:,:,first_period).^(1/theta_aux)))) / sum(sum(L_y_0));

pop_dyn(:,first_period) = pop_0;
pi_dyn(:,:,:,first_period) = pi_0;
u_dyn(:,first_period) = u_0;

for i=1:N
    for j=1:N
        prob_orig_given_dest_y_dyn(i,j,first_period) = sum(pi_dyn(i,j,:,first_period))*L_k_dyn(i,first_period) / sum(L_y_dyn(j,:,first_period));
    end
end
prob_orig_given_dest_y_dyn(:,:,first_period-1) = prob_orig_given_dest_y_dyn(:,:,first_period);

migration_exposure_dyn(:,:,first_period) = prob_orig_given_dest_y_dyn(:,:,first_period).*(1-eye(N)) ./ sum(prob_orig_given_dest_y_dyn(:,:,first_period).*(1-eye(N)),1);
migration_exposure_dyn(:,:,first_period-1) = migration_exposure_dyn(:,:,first_period);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% STEP 2: Dynamics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for t=first_period+1:first_period+T_aux-1

    t
    
    fertility_dyn(t) = (sum(data_pop(:,t)) - sum(L_k_dyn(:,t-1)) - sum(sum(L_y_dyn(:,:,t-1)))) / sum(L_k_dyn(:,t-1));
    adjustment_factor_0 = 1;
    
    data_pat_aux = data_pat(:,:,t);
    
    % Now, we do not know lambda_dyn yet. Assuming betta(BO) -> 1:
    % lambda_dyn = lambda_dyn_t-1 + SUM, where SUM = (Pat/Young)*G_t
    
    L_y_0 = (1/S)*(repmat(sum(L_y_dyn(:,:,t-1),2),1,S))*fertility_dyn(t); % Just a guess
    
    lambda_0 = lambda_dyn(:,:,t-1) + adjustment_factor_0*(data_pat_aux ./ L_y_0);
    
    u_0 = u_dyn(:,t-1);
    u_0 = u_0./geometric_mean(u_0);
    
    error_0 = 1000;
    
    while error_0 > 0.001
        
        kernel_utility_0 = u_0.*(lambda_0.^(1/theta_aux)); % The experience premium should appear as A_o^(1-betta) but does not matter, since it appears at the denominator too
        
        pi_0 = zeros(N,N,S);
    
        for j=1:N
            denom = sum(sum((kernel_utility_0./repmat(mu(j,:)',1,S)).^zeta));
            % Probability of going from j to any destination (:):
            pi_0(j,:,:) = ((kernel_utility_0./repmat(mu(j,:)',1,S)).^zeta)./denom;
        end
    
        L_y_1 = zeros(N,S);
        for j=1:N
            L_y_1 = L_y_1 + L_k_dyn(j,t-1)*squeeze(pi_0(j,:,:));
        end
        L_o_0 = L_y_dyn(:,:,t-1); L_k_0 = fertility_dyn(t)*sum(L_y_1,2);
        pop_0 = L_k_0 + sum(L_y_1,2) + sum(L_o_0,2);
        
        income_py_dyn(:,t) = sum(L_y_1.*(lambda_0.^(1/theta_aux)),2) ./ sum(L_y_1,2);
        aggr_income_py_dyn(t) = sum(sum(L_y_1.*(lambda_0.^(1/theta_aux)))) / sum(sum(L_y_1));
  
        aggregate_growth = aggr_income_py_dyn(1,t) / aggr_income_py_dyn(1,t-1);
        
        adjustment_factor_dyn(t) = (data_aggregate_growth / aggregate_growth)*adjustment_factor_0;
        
        u_1 = ((data_pop(:,t)./pop_0).^(1/zeta)).*u_0;
        u_1 = u_1./geometric_mean(u_1);
        
        error_0 = abs(adjustment_factor_0 - adjustment_factor_dyn(t)) + max(abs(u_0 - u_1)) + max(max(L_y_0 - L_y_1));
        
        rrr = rand;
        if rrr < 0.01
            error_0
        end

        adjustment_factor_0 = 0.99*adjustment_factor_0 + 0.01*adjustment_factor_dyn(t);
        
        L_y_0 = L_y_1;
        
        lambda_0 = zeros(N,S);
        for i=1:N
            for s=1:S
                if L_y_0(i,s) > 0
                    lambda_0(i,s) = delta*lambda_dyn(i,s,t-1) + adjustment_factor_dyn(t)*(data_pat_aux(i,s) / L_y_0(i,s));
                end
            end
        end
        
        u_0 = 0.99*u_0 + 0.01*u_1;
        
    end

    L_k_dyn(:,t) = L_k_0; L_y_dyn(:,:,t) = L_y_0; L_o_dyn(:,:,t) = L_o_0;
    pop_dyn(:,t) = pop_0;
    pi_dyn(:,:,:,t) = pi_0;
    u_dyn(:,t) = u_0;
    lambda_dyn(:,:,t) = lambda_0;
    income_pc_dyn(:,t) = (income_py_dyn(:,t).*sum(L_y_dyn(:,:,t),2) + A_o*income_py_dyn(:,t-1).*sum(L_o_dyn(:,:,t),2)) ./ (L_k_dyn(:,t) + sum(L_y_dyn(:,:,t),2) + sum(L_o_dyn(:,:,t),2));
    
    for i=1:N
        for j=1:N
            prob_orig_given_dest_y_dyn(i,j,t) = sum(pi_dyn(i,j,:,t))*L_k_dyn(i,t-1) / sum(L_y_dyn(j,:,t));
        end
    end
    migration_exposure_dyn(:,:,t) = prob_orig_given_dest_y_dyn(:,:,t).*(1-eye(N)) ./ sum(prob_orig_given_dest_y_dyn(:,:,t).*(1-eye(N)),1);

end

pop_growth_dyn = log(pop_dyn(:,2:end)) - log(pop_dyn(:,1:end-1));

% Weighted standard deviation in 1990 (in the data, 0.2144793
ave_log_income_1990 = (1/sum(pop_dyn(:,8)))*sum(log(income_pc_dyn(:,8)).*pop_dyn(:,8));
moments_vector(tthh) = sqrt((1/sum(pop_dyn(:,8)))*sum(((log(income_pc_dyn(:,8)) - ave_log_income_1990).^2).*pop_dyn(:,8)))

end

%%% Appendix Figure A.7:
plot(theta_vector,moments_vector,theta_vector,0.214*ones(1,length(theta_vector)),'k--')

csvwrite('results_model\show_identification.csv',moments_vector)
